#opções
options(scipen = 999)

#pacotes
library(haven)
library(psych)
library(memisc)
library(olsrr)
library(lm.beta)
library(mice)
library(e1071)
library(sjPlot)
library(sjlabelled)

#Importação das Bases de Dados

#Para a replicação deste script são necessários dois arquivos. O primeiro é o banco de dados
#integrado da sexta onda (versão 2018/09/12). Já o segundo é a base de dados longitudinal 
#do WVS (versão 2018/09/12). Ambos estão disponíveis em http://www.worldvaluessurvey.org/.

WV6 <- read_sav("WV6_Data_Spss_v20180912.sav")
alw6x <- subset(WV6, V2 == "32" | V2 == "76" | V2 == "152" | V2 == "170" |
                V2 == "484" | V2 == "604" | V2 == "858")
wvslong <- read_sav("WVS_Longitudinal_1981_2016_Spss_v20180912.sav")
wvslong <- subset(wvslong, S003 == "32" | S003 == "76" | S003 == "152" | S003 == "170" |
                    S003 == "484" | S003 == "604" | S003 == "858")
wvslong <- subset(wvslong, S002 == "5" | S002 == "6")
wvslong <- wvslong[c("S003", "S002", "E025", "E026", "E027", "E028", "E222")]
wvslong5 <- subset(wvslong, S002 == "5")
wvslong6 <- subset(wvslong, S002 == "6")

#Tabela 1 
#Os dados da Tabela 1 foram coletados a partir dos documentos técnicos
#das pesquisas realizados em cada país e se encontram em 
#http://www.worldvaluessurvey.org/WVSDocumentationWV6.jsp (Acesso em 15 de novembro de 2019).
#Para cada país coletaram-se os dados referentes ao tamanho amostral, margem de erro, 
#intervalo de confiança, nível de confiança e o ano em que o trabalho de camppo foi realizado.
#Todas as informações descritas acima foram agrupadas em uma Tabela, por país.

#Análise Fatorial
#Estala de protesto potencial (Tabela 2)
#ARG
tri.dep.ar <- alw6x[alw6x$V2 == "32", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.ar, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.ar)
alpha(tri.dep.ar)
#BRA
tri.dep.br <- alw6x[alw6x$V2 == "76", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.br, nfactors=1 ,n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.br)
alpha(tri.dep.br)
#CHI
tri.dep.ch <- alw6x[alw6x$V2 == "152", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.ch, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.ch)
alpha(tri.dep.ch)
#COL
tri.dep.co <- alw6x[alw6x$V2 == "170", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.co, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.co)
alpha(tri.dep.co)
#MEX
tri.dep.me <- alw6x[alw6x$V2 == "484", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.me, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.me)
alpha(tri.dep.me)
#PER
tri.dep.pe <- alw6x[alw6x$V2 == "604", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.pe, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.pe)
alpha(tri.dep.pe)
#URU
tri.dep.ur <- alw6x[alw6x$V2 == "858", c("V85", "V86", "V87", "V88", "V89")]
fa.poly(tri.dep.ur, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.dep.ur)
alpha(tri.dep.ur)

#Escala de Associativismo#
#ARG
tri.ass.ar <- alw6x[alw6x$V2 == "32", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.ar, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#BRA
tri.ass.br <- alw6x[alw6x$V2 == "76", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.br, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#CHI
tri.ass.ch <- alw6x[alw6x$V2 == "152", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.ch, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#COL
tri.ass.co <- alw6x[alw6x$V2 == "170", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.co, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#MEX
tri.ass.me <- alw6x[alw6x$V2 == "484", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.me, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#PER
tri.ass.pe <- alw6x[alw6x$V2 == "604", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.pe, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
#URU
tri.ass.ur <- alw6x[alw6x$V2 == "858", c("V26", "V27", "V28", "V29", "V30", "V31", "V32", "V33", "V34")]
fa.poly(tri.ass.ur, nfactors=3, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)

#Escala de Confiança Institucional (Tabela 3)#
#ARG
tri.cnf.ar <- alw6x[alw6x$V2 == "32", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.ar, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.ar)
alpha(tri.cnf.ar)
#BRA
tri.cnf.br <- alw6x[alw6x$V2 == "76", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.br, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.br)
alpha(tri.cnf.br)
#CHI
tri.cnf.ch <- alw6x[alw6x$V2 == "152", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.ch, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.ch)
alpha(tri.cnf.ch)
#COL
tri.cnf.co <- alw6x[alw6x$V2 == "170", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.co, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.co)
alpha(tri.cnf.co)
#MEX
tri.cnf.me <- alw6x[alw6x$V2 == "484", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.me, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.me)
alpha(tri.cnf.me)
#PER
tri.cnf.pe <- alw6x[alw6x$V2 == "604", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.pe, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.pe)
alpha(tri.cnf.pe)
#URU
tri.cnf.ur <- alw6x[alw6x$V2 == "858", c("V114", "V115", "V116", "V117")]
fa.poly(tri.cnf.ur, nfactors=1, n.obs = NA, rotate="oblimin", missing=TRUE, p =.05)
KMO(tri.cnf.ur)
alpha(tri.cnf.ur)

#Descritivas Medidas de Protesto (Tabela 4)
#5 onda
tab_xtab(wvslong5$S003, wvslong5$E025, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong5$S003, wvslong5$E026, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong5$S003, wvslong5$E027, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong5$S003, wvslong5$E028, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong5$S003, wvslong5$E222, show.row.prc = T, show.obs = F, show.na = T)
#6onda
tab_xtab(wvslong6$S003, wvslong6$E025, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong6$S003, wvslong6$E026, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong6$S003, wvslong6$E027, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong6$S003, wvslong6$E028, show.row.prc = T, show.obs = F, show.na = T)
tab_xtab(wvslong6$S003, wvslong6$E222, show.row.prc = T, show.obs = F, show.na = T)

#escalas - associativismo
alw6x$ass <- (alw6x$V26 + alw6x$V27 + alw6x$V28 + alw6x$V29 + alw6x$V30 + 
                alw6x$V31 + alw6x$V32 + alw6x$V33 + alw6x$V34)
alw6x$ass <- as.numeric(alw6x$ass)
alw6x$ass <- recode(alw6x$ass, 0 <- 0, 1 <- c(1:18))

#escalas - confiança institucional
alw6x$V114 <- as.factor(alw6x$V114)
alw6x$V115 <- as.factor(alw6x$V115)
alw6x$V116 <- as.factor(alw6x$V116)
alw6x$V117 <- as.factor(alw6x$V117)
alw6x$V114 <- recode(alw6x$V114, 1 <- 4, 2 <- 3, 3 <- 2, 4 <- 1)
alw6x$V115 <- recode(alw6x$V115, 1 <- 4, 2 <- 3, 3 <- 2, 4 <- 1)
alw6x$V116 <- recode(alw6x$V116, 1 <- 4, 2 <- 3, 3 <- 2, 4 <- 1)
alw6x$V117 <- recode(alw6x$V117, 1 <- 4, 2 <- 3, 3 <- 2, 4 <- 1)
alw6x$V114 <- as.numeric(alw6x$V114)
alw6x$V115 <- as.numeric(alw6x$V115)
alw6x$V116 <- as.numeric(alw6x$V116)
alw6x$V117 <- as.numeric(alw6x$V117)
alw6x$conf.inx <- alw6x$V114 + alw6x$V115 + alw6x$V116 + alw6x$V117
alw6x$conf.in <- (alw6x$conf.inx - 4)
summary(alw6x$conf.in)

#escalas - protesto potencial
alw6x$V85 <- as.factor(alw6x$V85)
alw6x$V86 <- as.factor(alw6x$V86)
alw6x$V87 <- as.factor(alw6x$V87)
alw6x$V88 <- as.factor(alw6x$V88)
alw6x$V89 <- as.factor(alw6x$V89)
alw6x$V85 <- recode(alw6x$V85, 1 <- 3, 2 <- 2, 3 <- 1)
alw6x$V86 <- recode(alw6x$V86, 1 <- 3, 2 <- 2, 3 <- 1)
alw6x$V87 <- recode(alw6x$V87, 1 <- 3, 2 <- 2, 3 <- 1)
alw6x$V88 <- recode(alw6x$V88, 1 <- 3, 2 <- 2, 3 <- 1)
alw6x$V89 <- recode(alw6x$V89, 1 <- 3, 2 <- 2, 3 <- 1)
alw6x$V85 <- as.numeric(alw6x$V85)
alw6x$V86 <- as.numeric(alw6x$V86)
alw6x$V87 <- as.numeric(alw6x$V87)
alw6x$V88 <- as.numeric(alw6x$V88)
alw6x$V89 <- as.numeric(alw6x$V89)
alw6x$protpon <- ((alw6x$V85 + alw6x$V86 + alw6x$V87 + alw6x$V88 + alw6x$V89)-5)
summary(alw6x$protpon)
alw6x$protpon <- as.numeric(alw6x$protpon)

#recodificação de variável - escolaridade
alw6x$V248 <- as.factor(alw6x$V248)
alw6x$esc<- recode(alw6x$V248, 1 <- c("1", "2", "3", "4"), 2 <- c("5", "6", "7"), 
                   3 <- c("8", "9"))

#recodificação de variável - sexo
alw6x$sexo <- alw6x$V240

#recodificação de variável - idade
alw6x$idade <- alw6x$V242

#reodificação de variável - renda
alw6x$renda <- alw6x$V239

#recodificação de variável - Valores Emancipatórios
alw6x$EV <- alw6x$RESEMAVAL

#redução BD
BD <- alw6x[c("V2", "ass", "conf.in", "renda", "esc", "sexo",
              "idade", "EV", "protpon", "V23")]

#imputação de casos faltantes
BDimp <- mice(BD,m=5,maxit=50,meth='pmm',seed=500)
BDimp2 <- complete(BDimp, 1)

#redução BD - Países
#argentina
BDar <- BDimp2[BDimp2$V2 == "32", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                    "idade", "EV", "protpon", "V23")]
#Brasil
BDbr <- BDimp2[BDimp2$V2 == "76", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                    "idade", "EV", "protpon", "V23")]
#Chile
BDch <- BDimp2[BDimp2$V2 == "152", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                     "idade", "EV", "protpon", "V23")]

#Colômbia
BDco <- BDimp2[BDimp2$V2 == "170", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                     "idade", "EV", "protpon", "V23")]

#Mexico
BDme <- BDimp2[BDimp2$V2 == "484", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                     "idade", "EV", "protpon", "V23")]

#Peru
BDpe <- BDimp2[BDimp2$V2 == "604", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                     "idade", "EV", "protpon", "V23")]

#Uruguai
BDur <- BDimp2[BDimp2$V2 == "858", c("V2", "ass", "conf.in", "renda", "esc", "sexo",
                                     "idade", "EV", "protpon", "V23")]

#Modelos Lineares (Tabela 5)
#modelo linear geral (todos os casos)
BDimp2$V2 <- as.factor(BDimp2$V2)
BDimp2$V23 <- as.numeric(BDimp2$V23)
modelo1 <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                EV + V2, data = BDimp2)

summary(modelo1)
ols_vif_tol(modelo1)
ols_eigen_cindex(modelo1)
modelo1.beta <- lm.beta(modelo1)
summary(modelo1.beta)
plot(modelo1)


#modelo linear Argentina
modelo1ar <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDar)

summary(modelo1ar)
ols_vif_tol(modelo1ar)
ols_eigen_cindex(modelo1ar)
modelo1ar.beta <- lm.beta(modelo1ar)
summary(modelo1ar.beta)
plot(modelo1ar)


#modelo linear brasil
modelo1br <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDbr)

summbry(modelo1br)
ols_vif_tol(modelo1br)
ols_eigen_cindex(modelo1br)
modelo1br.beta <- lm.beta(modelo1br)
summary(modelo1br.beta)
plot(modelo1br)


#modelo linear Chile
modelo1ch <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDch)

summary(modelo1ch)
ols_vif_tol(modelo1ch)
ols_eigen_cindex(modelo1ch)
modelo1ch.beta <- lm.beta(modelo1ch)
summary(modelo1ch.beta)
plot(modelo1ch)


#modelo linear Colombia
modelo1co <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDco)

summary(modelo1co)
ols_vif_tol(modelo1co)
ols_eigen_cindex(modelo1co)
modelo1co.beta <- lm.beta(modelo1co)
summary(modelo1co.beta)
plot(modelo1co)


#modelo linear México
modelo1me <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDme)

summmey(modelo1me)
ols_vif_tol(modelo1me)
ols_eigen_cindex(modelo1me)
modelo1me.beta <- lm.beta(modelo1me)
summary(modelo1me.beta)
plot(modelo1me)


#modelo linear Peru
modelo1pe <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDpe)

summary(modelo1pe)
ols_vif_tol(modelo1pe)
ols_eigen_cindex(modelo1pe)
modelo1pe.beta <- lm.beta(modelo1pe)
summary(modelo1pe.beta)
plot(modelo1pe)

#modelo linear Uruguai
modelo1ur <- lm(protpon ~ sexo + idade + esc + renda + ass + V23 + conf.in  + 
                  EV, data = BDur)

summury(modelo1ur)
ols_vif_tol(modelo1ur)
ols_eigen_cindex(modelo1ur)
modelo1ur.beta <- lm.beta(modelo1ur)
summary(modelo1ur.beta)
plot(modelo1ur)

#Betas padronizados (Tabela 6)

modelo1.beta <- lm.beta(modelo1)
modelo1ar.beta <- lm.beta(modelo1ar)
modelo1br.beta <- lm.beta(modelo1br)
modelo1ch.beta <- lm.beta(modelo1ch)
modelo1co.beta <- lm.beta(modelo1co)
modelo1me.beta <- lm.beta(modelo1me)
modelo1pe.beta <- lm.beta(modelo1pe)
modelo1ur.beta <- lm.beta(modelo1ur)

#Tabulação dos Modelos#
#Tabela 5
tab_model(modelo1, modelo1ar, modelo1br, modelo1ch, modelo1co, modelo1me,
          modelo1pe, modelo1ur, show.ci = F, collapse.se = T, p.style = "asterisk")

#Tabela 6 (Beta Padronizado)
tab_model(modelo1.beta, modelo1ar.beta, modelo1br.beta, modelo1ch.beta, modelo1co.beta,
          modelo1me.beta, modelo1pe.beta, modelo1ur.beta, 
          show.ci = F, collapse.se = T, p.style = "asterisk")

#Gráficos
#O Gráfico 1 utiliza dados provenientes do World Values Survey e do Varieties of Democracy.
#O primeiro já foi importado acima e pode ser acessado em http://www.worldvaluessurvey.org/
#(versão 2018/09/12). Já o segundo está disponpivel em 
#https://www.v-dem.net/en/data/data-version-10/, (versão 10). 

vdem <- read_sav("V-Dem-CY-Core-v10.sav")
vdema <- subset(vdem, vdem$country_name == "Argentina" | vdem$country_name == "Brazil" |
                vdem$country_name == "Chile" | vdem$country_name == "Colombia" |
                vdem$country_name == "Mexico" | vdem$country_name == "Peru" | 
                vdem$country_name == "Uruguay")
vdema <- subset(vdema, vdema$year == "2013")
vdema <- vdema[c("country_name", "v2x_libdem")]
vdema <- vdema[order(vdema$country_name),]
BDimp2$protpong <- BDimp2$protpon / 10
ev <- aggregate(BDimp2["EV"], BDimp2["V2"], FUN = mean)
pp <- aggregate(BDimp2["protpong"], BDimp2["V2"], FUN = mean)
graf1 <- cbind(vdema, ev, pp)
write.csv2(graf1, "graf1.csv")

#Depois de salvo, o arquivo graf1.cvs foi importado para o Excel. A variável "v2x_libdem"
#foi renomeada para "Qualidade da Democracia", a variável "EV" passou a se chamar 
#"Valores Emancipatório" e, por sua vez, va variável "protpong" foi renomeada para 
#"Protesto Potencial". O Gráfico 1 foi plotado no Excel.

#Análise de Resíduos x Valores Previstos (Gráfico 2)
windowsFonts(A = windowsFont("Cambria"))
par(mfrow=c(2,4))
plot(x=fitted(modelo1), y=residuals(modelo1), main = "All Countries", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1ar), y=residuals(modelo1ar), main = "Argentina", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1br), y=residuals(modelo1br), main = "Brazil", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1ch), y=residuals(modelo1ch), main = "Chile",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1co), y=residuals(modelo1co), main = "Colombia", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1me), y=residuals(modelo1me), main = "Mexico", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1pe), y=residuals(modelo1pe), main = "Peru", family = "A",
     xlab = "fitted", ylab = "residuals")
plot(x=fitted(modelo1ur), y=residuals(modelo1ur), main = "Uruguai", family = "A",
     xlab = "fitted", ylab = "residuals")


#QQ-plot (Gráfico 3)
par(mfrow=c(2,4))
qqnorm(rstandard(modelo1), main = "All Countries", family = "A")
qqline(rstandard(modelo1))
qqnorm(rstandard(modelo1ar), main = "Argentina", family = "A")
qqline(rstandard(modelo1ar))
qqnorm(rstandard(modelo1br), main = "Brazil", family = "A")
qqline(rstandard(modelo1br))
qqnorm(rstandard(modelo1ch), main = "Chile", family = "A")
qqline(rstandard(modelo1ch))
qqnorm(rstandard(modelo1co), main = "Colombia", family = "A")
qqline(rstandard(modelo1co))
qqnorm(rstandard(modelo1me), main = "Mexico", family = "A")
qqline(rstandard(modelo1me))
qqnorm(rstandard(modelo1pe), main = "Peru", family = "A")
qqline(rstandard(modelo1pe))
qqnorm(rstandard(modelo1ur), main = "Uruguai", family = "A")
qqline(rstandard(modelo1ur))

##FIM DO SCRIPT##

